function penloglik = cicaar_bayes(D,KC,L,M,interceptflag,ll,Hinv)
% CICAAR_BAYES   evaluate model posterior of CICAAR decomposition
%
%    penloglik = cicaar_bayes(D,KC,L,M,interceptflag,ll,Hinv) % Laplace
%    penloglik = cicaar_bayes(D,KC,L,M,interceptflag,ll,N)    % BIC
%
% inputs:
%
%    D,KC,L,M  :  CICAAR parameters
%  interceptflag  :  do. 0: no intercept, 1: intercept
%             ll  :  returned loglik from CICAAR
%           Hinv  :  inverse Hessian at optimum in PACKED TRIANGULAR
%              N  :  number of samples
%
% output:
%
%    penloglik  :  penalized log likelihood (the higher the better)
%
% Copyright: Mads Dyrholm, 2007
interceptflag = (interceptflag~=0);
dim = D*D + D*KC*L + D*M + interceptflag*D;
if size(Hinv) == [1 1]
  % BIC
  bic = - dim/2 * log(Hinv);
  fprintf('BIC term added: %f\n',bic);
  penloglik = ll + bic;
else
  % Laplace
  Hinvfull = ones(dim,dim);
  Hinvfull = tril(Hinvfull);
  Hinvfull(find(Hinvfull)) = Hinv;
  Hinvfull = Hinvfull + tril(Hinvfull,-1)'; % full hessian
  U = chol(Hinvfull);             %  stable method for log determinant
  logdet = 2*sum(log(diag(U)));   %  of SPD matrix
  laplace = dim/2*log(2*pi) + 0.5*logdet;
  fprintf('Laplace approcimation used, added: %g\n',laplace);
  penloglik = ll + laplace;
end

  
 
